###############################################################################################################################################
### Setup

# Set working directory
setwd('~/Dropbox/Projects/Minority Party/Paper 2/APSR Replication Files/')

# Source setup file
source('setup.R')


### Load bills data
bills <- fread(file='bills97-109_partyPriorities_issuePreferences.csv', data.table=FALSE, stringsAsFactors = FALSE) %>%
  filter(congress %in% 99:109)


### Add in SOTU data
# Load, filter to years matching the congresses studied (1985-2006, though the 109th technically ends in 2007 it was before the 2007 sotu)
sotu <- fread(file='SOTU-topicmentions-1946-2019.csv', data.table=FALSE, stringsAsFactors=FALSE) %>%
  filter(year %in% 1985:2006) %>%
  dplyr::rename(minor=subtopic) 

# Match years to congresses
sotu$congress <- NA
congs <- sort(unique(bills$congress))
years <- sort(unique(sotu$year))
for(ii in 1:d(years)){
  if(years[ii] %% 2 == 1){ # If it's an odd year
    sotu$congress[sotu$year==years[ii]] <- congs[(ii+1)/2]
  } else {
    sotu$congress[sotu$year==years[ii]] <- congs[(ii/2)]
  }
}

# Drop years
sotu %<>% dplyr::select(-year)

# Join with bills
bills <- left_join(bills, sotu, by=c('congress', 'minor'))


### Miscellaneous formatting

# Add in 0s for sotu variables (issues that weren't mentioned)
bills$any_statement[is.na(bills$any_statement) & !is.na(bills$minor)] <- 0
bills$num_statements[is.na(bills$num_statements) & !is.na(bills$minor)] <- 0
bills$perc_statements[is.na(bills$perc_statements) & !is.na(bills$minor)] <- 0
bills$perc_topic_statements[is.na(bills$perc_topic_statements) & !is.na(bills$minor)] <- 0

# Cut down to distinct rows (if any duplicates)
bills %<>% distinct()

# Indicator for divided government
bills %<>% mutate(divided=ifelse(!congress %in% c(103, 107, 108, 109), 1, 0))

# Remake majPriorityAll and minPriorityAll into factors (makes plotting better)
bills$majPriorityAll <- factor(bills$majPriorityAll, labels=c('No', 'Yes')) %>% relevel(., ref='No')
bills$minPriorityAll <- factor(bills$minPriorityAll, labels=c('No', 'Yes')) %>% relevel(., ref='No')
bills$any_statement <- factor(bills$any_statement, labels=c('No', 'Yes')) %>% relevel(., ref='No')

# Filter to substantive bills in the 99th-109th congresses (dataframe for modeling)
mdat <- filter(bills, (congress %in% 99:109) & substantive==1)



###############################################################################################################################################
### Figure A1: Correlation between covariates

### Set up data

# Vectors of IVs and covariates
covs <- c('leader', 'les', 'mref', 'absPVI', 
          'closerMinCS', 'elec.year', 'relSen', 'medDiff', 'majPriorityAll',
          'sizeMajH', 'numSenate', 'unityPty', 'presControl') 
mainIV <- c('majSD', 'minSD', 'any_statement') 


### Check for multicollinearity in predictors

# Pull correlations between all IVs
ivMat <- dplyr::select(mdat, all_of(c(covs, mainIV))) %>%
  mutate(any_statement=car::recode(any_statement, "'Yes'=1; 'No'=0") %>% num()) %>%
  mutate(majPriorityAll=car::recode(majPriorityAll, "'Yes'=1; 'No'=0") %>% num())
corMat <- cor(ivMat, use='pairwise.complete.obs')

# Format names for plotting
rownames(corMat) <- colnames(corMat) <- c('Leader', 'Legislative Effectiveness', 'Member Referral Committee', 'Electoral Safety', 'Closer to Minority',
                                          'Election Year', 'Related in Senate', 'Party Difference', 'Majority Priority', 'Size House Majority', 
                                          'Seats in Senate', 'Party Unity', 'WH Control', 'Majority Spread', 'Minority Spread', 'Minority Priority')

# Plot correlations between IVs
pdf(file='figureA1.pdf', width=10, height=10)
corrplot::corrplot(corMat, method='circle', tl.col='black', type='upper')
dev.off()



############################################################################################################################################################################
### Table A1: Models building from main IVs to full model

# Just with the operationalizations of main IVs
m.theory <- glm(pv_all ~ majSD * minSD * any_statement, 
                family=binomial, data=filter(mdat, divided==1 & congress %in% 99:103))

# Add sponsor controls
m.spon <- glm(pv_all ~ majSD * minSD * any_statement +
                majority + leader + les + mref + absPVI, 
              family=binomial, data=filter(mdat, divided==1 & congress %in% 99:103))

# Add bill controls
m.bill <- glm(pv_all ~ majSD * minSD * any_statement +
                majority + leader + les + mref + absPVI + 
                closerMinCS + elec.year + relSen + medDiff + majPriorityAll, 
              family=binomial, data=filter(mdat, divided==1 & congress %in% 99:103))

# Full model
m.full <- glm(pv_all ~ majSD * minSD * any_statement + 
                majority + leader + les + mref + absPVI + 
                closerMinCS + elec.year + relSen + medDiff + majPriorityAll + 
                numSenate + sizeMajH, 
              family=binomial, data=filter(mdat, divided==1 & congress %in% 99:103))

# R-squareds for table
r2s.build <- sapply(list(m.theory, m.spon, m.bill, m.full), 
                    function(x){
                      round(NagelkerkeR2(x)$R2, 3)
                    })

# Save model output
stargazer(m.theory, m.spon, m.bill, m.full, type='latex', out='tableA1.tex',
          column.sep.width='0pt', no.space=TRUE, font.size='tiny', dep.var.labels.include = FALSE, dep.var.caption='', 
          title='Logistic regression models predicting whether a bill receives a final passage vote in the House, building from main IVs to the full model.', model.numbers=FALSE,
          column.labels = c('Main IVs', 'Sponsor Controls', 'Bill Controls', 'Full Model'),
          covariate.labels = c( # Main IV labels
            'Majority Spread', 'Minority Spread', 'Minority Priority',
            # Sponsor covariate labels
            'Majority', 'Leadership', 'Legislative Effectiveness', 'Member Referral Committee', 'Electoral Safety', 
            # Bill covariate labels
            'Closer to Minority', 'Election Year', 'Related in Senate', 'Party Difference', 'Majority Priority',
            # Political environment covariate labels
            'Seats in Senate', 'House Majority Size',
            # Interaction labels
            'Majority Spread * Minority Spread', 'Majority Spread * SOTU Mention', 'Minority Spread * SOTU Mention', 
            'Majority Spread * Minority Spread * SOTU Mention'
          ),
          keep.stat='n', add.lines = list(c('Nagelkerke Pseudo R2', r2s.build))
)



################################################################################################################################################################
### Figures A2-A4: Fixed effects for Congress

### Run models

# Full model (opportunity and cohesion), no interaction
m.noint.fe <- glm(pv_all ~ majSD + minSD + 
                 majority + leader + les + mref + absPVI + 
                 closerMinCS + elec.year + relSen + majPriorityAll +
                 factor(congress) - 1, 
               family=binomial, data=mdat)

# Full model (add interaction)
m.full.fe <- glm(pv_all ~ majSD * minSD + 
                majority + leader + les + mref + absPVI + 
                closerMinCS + elec.year + relSen + majPriorityAll + 
                factor(congress) - 1, 
              family=binomial, data=mdat)

# Model including SOTU mentions--adds motivation analysis 
# (only in cases of divided government between House and White House, which cuts out 103rd, 107th-109th congresses)
m.sotu.fe <- glm(pv_all ~ majSD * minSD * any_statement + 
                majority + leader + les + mref + absPVI + 
                closerMinCS + elec.year + relSen + majPriorityAll +
                factor(congress) - 1, 
              family=binomial, 
              data=filter(mdat, divided==1))

# Models computing whether a bill became law (among all bills and those that passed the House)
covs <- ' + majority + leader + les + mref + absPVI +
                closerMinCS + elec.year + relSen + majPriorityAll +
                factor(congress) - 1'
m.law.fe <- glm(paste0('plaw ~ majSD * minSD * any_statement', covs),
             family=binomial, data=filter(mdat, divided==1))
m.law.pass.fe <- glm(paste0('plaw ~ majSD * minSD * any_statement', covs),
                  family=binomial, data=filter(mdat, divided==1 & passh==1))

# Define plot colors
cols3 <- brewer.pal(6, 'Greys') %>% .[c(4:6)]



### Figure A2: 2-way interaction (majority spread and minority spread)

# Get data
pDat <- predDat(model=m.full.fe, vars=c('minSD [all]', 'majSD'), data=mdat, std=TRUE)

# Plot
pdf("figureA2.pdf", width=5, height=5)

breaks <- 3
plot(0:(nrow(pDat)+2*breaks-1), 0:(nrow(pDat)+2*breaks-1), type='n', 
     ylab='Probability of Vote', xlab='Minority Spread', 
     ylim=c(0.02,0.07), xaxt='n', las=2, main='')
for(ii in 1:nrow(pDat)){
  colIndex <- which(c('Low', 'Median', 'High')==pDat$majSD[ii])
  if(ii %in% 1:3){
    points(ii, pDat$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii, 2), y=c(pDat$lo[ii], pDat$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 4:6){
    points(ii+2, pDat$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+2, 2), y=c(pDat$lo[ii], pDat$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 7:9) {
    points(ii+4, pDat$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+4, 2), y=c(pDat$lo[ii], pDat$hi[ii]), col=cols3[colIndex], lwd=3)
  }
}
axis(side=1, at=c(2, 7, 12), labels=c('Low', 'Median', 'High'), las=1, tick=FALSE)
legend(x='topright', legend=c('Low', 'Median', 'High'), fill=cols3, bty='n', title='Majority Spread')

dev.off()



### Figure A3: 3-way interaction

# Get data
p3 <- predDat(model=m.sotu.fe, vars=c('any_statement', 'minSD', 'majSD'))

# Set up plot parameters
pdf("figureA3.pdf", width=10, height=5)
par(mfrow=c(1,2))

p3No <- filter(p3, x==1)
breaks <- 3
plot(0:(nrow(pDat)+2*breaks-1), 0:(nrow(pDat)+2*breaks-1), type='n', 
     ylab='Probability of Vote', xlab='Minority Spread', 
     ylim=c(0.01, 0.06), xaxt='n', las=2, main='Not a Minority Priority')
for(ii in 1:nrow(p3No)){
  colIndex <- which(c('Low', 'Median', 'High')==p3No$majSD[ii])
  if(ii %in% 1:3){
    points(ii, p3No$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii, 2), y=c(p3No$lo[ii], p3No$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 4:6){
    points(ii+2, p3No$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+2, 2), y=c(p3No$lo[ii], p3No$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 7:9) {
    points(ii+4, p3No$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+4, 2), y=c(p3No$lo[ii], p3No$hi[ii]), col=cols3[colIndex], lwd=3)
  }
}
axis(side=1, at=c(2, 7, 12), labels=c('Low', 'Median', 'High'), las=1, tick=FALSE)
legend(x='topright', legend=c('Low', 'Median', 'High'), fill=cols3, bty='n', title='Majority Spread')

p3Yes <- filter(p3, x==2)
breaks <- 3
plot(0:(nrow(pDat)+2*breaks-1), 0:(nrow(pDat)+2*breaks-1), type='n', 
     ylab='Probability of Vote', xlab='Minority Spread', 
     ylim=c(0.01,0.06), xaxt='n', las=2, main='Minority Priority')
for(ii in 1:nrow(p3Yes)){
  colIndex <- which(c('Low', 'Median', 'High')==p3Yes$majSD[ii])
  if(ii %in% 1:3){
    points(ii, p3Yes$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii, 2), y=c(p3Yes$lo[ii], p3Yes$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 4:6){
    points(ii+2, p3Yes$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+2, 2), y=c(p3Yes$lo[ii], p3Yes$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 7:9) {
    points(ii+4, p3Yes$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+4, 2), y=c(p3Yes$lo[ii], p3Yes$hi[ii]), col=cols3[colIndex], lwd=3)
  }
}
axis(side=1, at=c(2, 7, 12), labels=c('Low', 'Median', 'High'), las=1, tick=FALSE)
legend(x='topleft', legend=c('Low', 'Median', 'High'), fill=cols3, bty='n', title='Majority Spread')

par(mfrow=c(1,1))
dev.off()


### Figure A4: Probability of becoming law
pLaw <- predDat(m.law.fe, vars=c('any_statement', 'minSD', 'majSD'))

pdf("figureA4.pdf", width=10, height=5)
par(mfrow=c(1,2))

pLawNo <- filter(pLaw, x==1)
breaks <- 3
plot(0:(nrow(pLawNo)+2*breaks), 0:(nrow(pLawNo)+(2*breaks)), type='n', 
     ylab='Probability of Becoming Law', xlab='Minority Spread', 
     ylim=c(0,0.04), xaxt='n', las=2, main='Not a Minority Priority')
for(ii in 1:nrow(pLawNo)){
  colIndex <- which(c('Low', 'Median', 'High')==pLawNo$majSD[ii])
  if(ii %in% 1:3){
    points(ii, pLawNo$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii, 2), y=c(pLawNo$lo[ii], pLawNo$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 4:6){
    points(ii+2, pLawNo$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+2, 2), y=c(pLawNo$lo[ii], pLawNo$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 7:9) {
    points(ii+4, pLawNo$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+4, 2), y=c(pLawNo$lo[ii], pLawNo$hi[ii]), col=cols3[colIndex], lwd=3)
  }
}
axis(side=1, at=c(2, 7, 12), labels=c('Low', 'Median', 'High'), las=1, tick=FALSE)
legend(x='topright', legend=c('Low', 'Median', 'High'), fill=cols3, bty='n', title='Majority Spread')

pLawYes <- filter(pLaw, x==2)
breaks <- 3
plot(0:(nrow(pLawYes)+2*breaks), 0:(nrow(pLawYes)+(2*breaks)), type='n', 
     ylab='Probability of Becoming Law', xlab='Minority Spread', 
     ylim=c(0,0.04), xaxt='n', las=2, main='Minority Priority')
for(ii in 1:nrow(pLawYes)){
  colIndex <- which(c('Low', 'Median', 'High')==pLawYes$majSD[ii])
  if(ii %in% 1:3){
    points(ii, pLawYes$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii, 2), y=c(pLawYes$lo[ii], pLawYes$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 4:6){
    points(ii+2, pLawYes$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+2, 2), y=c(pLawYes$lo[ii], pLawYes$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 7:9) {
    points(ii+4, pLawYes$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+4, 2), y=c(pLawYes$lo[ii], pLawYes$hi[ii]), col=cols3[colIndex], lwd=3)
  }
}
axis(side=1, at=c(2, 7, 12), labels=c('Low', 'Median', 'High'), las=1, tick=FALSE)
legend(x='topright', legend=c('Low', 'Median', 'High'), fill=cols3, bty='n', title='Majority Spread')

par(mfrow=c(1,1))
dev.off()



#########################################################################################################################################################################
### Table A2 and Figure A5: Bills by party of introducing member

### Set up data

# Vectors of IVs and covariates
# Vector of covariates for formulae
covs <- ' + majority + leader + les + mref + absPVI +
closerMinCS + elec.year + relSen + medDiff + majPriorityAll +
numSenate + sizeMajH'

# Model including SOTU mentions--adds motivation analysis 
m.sotu <- glm(paste0('pv_all ~ majSD * minSD * any_statement', covs), 
              family=binomial, data=filter(mdat, divided==1))
m.sotu.min <- glm(paste0('pv_all ~ majSD * minSD * any_statement', covs), 
                  family=binomial, data=filter(mdat, divided==1 & majority==0))
m.sotu.maj <- glm(paste0('pv_all ~ majSD * minSD * any_statement', covs), 
                  family=binomial, data=filter(mdat, divided==1 & majority==1))

# Model Table 
stargazer(m.sotu, m.sotu.min, m.sotu.maj, type='latex', out='tableA2.tex',
          column.sep.width='0pt', no.space=TRUE, font.size='tiny', dep.var.labels.include = FALSE, dep.var.caption='', 
          title='Logistic regression models predicting whether a bill receives a final passage vote in the House by party of the introducing member.', model.numbers=FALSE,
          column.labels = c('All Bills', 'Minority Sponsor', 'Majority Sponsor'),
          covariate.labels = c( # Main IV labels
            'Majority Spread', 'Minority Spread', 'Minority Priority',
            # Sponsor covariate labels
            'Majority', 'Leadership', 'Legislative Effectiveness', 'Member Referral Committee', 'Electoral Safety', 
            # Bill covariate labels
            'Closer to Minority', 'Election Year', 'Related in Senate', 'Party Difference', 'Majority Priority',
            # Political environment covariate labels
            'Seats in Senate', 'House Majority Size',
            # Interaction labels
            'Majority Spread * Minority Spread', 'Majority Spread * Minority Priority', 'Minority Spread * Minority Priority', 
            'Majority Spread * Minority Spread * Minority Priority'
          ),
          keep.stat='n'
)

### Figure A5: 3-way interaction for all bills and minority sponsored bills

# Data
pMin <- predDat(model=m.sotu.min, vars=c('any_statement', 'minSD', 'majSD'))
pSotu <- predDat(model=m.sotu, vars=c('any_statement', 'minSD', 'majSD'))

# Plot
pdf("figureA5.pdf", width=10, height=10)
par(mfrow=c(2,2))
predProbPlot(pDat = pSotu, breaks=3, ylab = 'Probability of Vote', xlab = 'Minority Spread', ylim=c(0, 0.06), var='SD',
             mains=c('Not a Minority Priority', 'Minority Priority'), rowControl = NULL, legTitle = 'Majority Spread')
predProbPlot(pDat = pMin, breaks = 3, ylab = 'Probability of Vote', xlab = 'Minority Spread', ylim=c(0, 0.06), var='SD',
             mains=c('Not a Minority Priority (Minority Sponsor)', 'Minority Priority (Minority Sponsor)'), rowControl = NULL, legTitle = 'Majority Spread')
par(mfrow=c(1,1))
dev.off()



################################################################################################################################################################
### Figures A6-A8: Models separately by period

### Covariates
covs <- ' + majority + leader + les + mref + absPVI + closerMinCS + elec.year + 
relSen + medDiff + majPriorityAll + numSenate + sizeMajH'


### Run models

# Democratic majority (99th-103rd)
m.sotu.d <- glm(paste0('pv_all ~ majSD * minSD * any_statement + ', covs), 
                family=binomial, data=filter(mdat, divided==1 & congress %in% 99:103))
m.law.d <- glm(paste0('plaw ~ majSD * minSD * any_statement +', covs),
               family=binomial, data=filter(mdat, divided==1 & congress %in% 99:103))
m.law.pass.d <- glm(paste0('plaw ~ majSD * minSD * any_statement', covs),
                    family=binomial, data=filter(mdat, divided==1 & passh==1 & congress %in% 99:103))


# Republican majority (94th-109th)
m.sotu.r <- glm(paste0('pv_all ~ majSD * minSD * any_statement + ', covs), 
                family=binomial, data=filter(mdat, divided==1 & congress %in% 104:109))
m.law.r <- glm(paste0('plaw ~ majSD * minSD * any_statement +', covs),
               family=binomial, data=filter(mdat, divided==1 & congress %in% 104:109))
m.law.pass.r <- glm(paste0('plaw ~ majSD * minSD * any_statement', covs),
                    family=binomial, data=filter(mdat, divided==1 & passh==1 & congress %in% 104:109))


### Figures

# Get plot data
pSotuD <- predDat(model=m.sotu.d, vars=c('any_statement', 'minSD', 'majSD')) 
pSotuR <- predDat(model=m.sotu.r, vars=c('any_statement', 'minSD', 'majSD')) 
pLawD <- predDat(model=m.law.d, vars=c('any_statement', 'minSD', 'majSD')) 
pLawR <- predDat(model=m.law.r, vars=c('any_statement', 'minSD', 'majSD')) 
pLawPassD <- predDat(model=m.law.pass.d, vars=c('any_statement', 'minSD', 'majSD')) 
pLawPassR <- predDat(model=m.law.pass.r, vars=c('any_statement', 'minSD', 'majSD')) 


### Figure A6: 3-way interaction, Pr(Vote), by period
# Democratic majority
pdf("figureA6_dem.pdf", width=10, height=5)
predProbPlot(pDat = pSotuD, breaks = 3, ylab = 'Probability of Vote', xlab = 'Minority Spread', ylim=c(0.01, 0.08), var='SD',
             mains=c('Not a Minority Priority', 'Minority Priority'), rowControl = c(1,2), legTitle = 'Majority Spread')
dev.off()

# Republican majority
pdf("figureA6_rep.pdf", width=10, height=5)
predProbPlot(pDat = pSotuR, breaks = 3, ylab = 'Probability of Vote', xlab = 'Minority Spread', ylim=c(0.01, 0.08), var='SD',
             mains=c('Not a Minority Priority', 'Minority Priority'), rowControl = c(1,2), legTitle = 'Majority Spread')
dev.off()


### Figure A7: 3-way interaction, Pr(law), by period
# Democratic majority
pdf("figureA7_dem.pdf", width=10, height=5)
predProbPlot(pDat = pLawD, breaks = 3, ylab = 'Probability of Becoming Law', xlab = 'Minority Spread', ylim=c(0, 0.05), var='SD',
             mains=c('Not a Minority Priority', 'Minority Priority'), rowControl = c(1,2), legTitle = 'Majority Spread')
dev.off()

# Republican majority
pdf("figureA7_rep.pdf", width=10, height=5)
predProbPlot(pDat = pLawR, breaks = 3, ylab = 'Probability of Becoming Law', xlab = 'Minority Spread', ylim=c(0, 0.05), var='SD',
             mains=c('Not a Minority Priority', 'Minority Priority'), rowControl = c(1,2), legTitle = 'Majority Spread')
dev.off()


### Figure A8: 3-way interaction, Pr(law|passed house), by period
# Democratic majority
pdf("figureA8_dem.pdf", width=10, height=5)
predProbPlot(pDat = pLawPassD, breaks = 3, ylab = 'Probability of Becoming Law', xlab = 'Minority Spread', ylim=c(0.25, 0.65), var='SD',
             mains=c('Not a Minority Priority', 'Minority Priority'), rowControl = c(1,2), legTitle = 'Majority Spread')
dev.off()

# Republican majority
pdf("figureA8_rep.pdf", width=10, height=5)
predProbPlot(pDat = pLawPassR, breaks = 3, ylab = 'Probability of Becoming Law', xlab = 'Minority Spread', ylim=c(0.25, 0.65), var='SD',
             mains=c('Not a Minority Priority', 'Minority Priority'), rowControl = c(1,2), legTitle = 'Majority Spread')
dev.off()



###############################################################################################################################################################
### Table A3 and Figure A9: Cutting public lands bills

### Run models

# Full models
m.full <- glm(pv_all ~ majSD * minSD * any_statement + 
                majority + leader + les + mref + absPVI + 
                closerMinCS + elec.year + relSen + medDiff + majPriorityAll + 
                numSenate + sizeMajH, 
              family=binomial, data=filter(mdat, divided==1))


# Model cutting out public lands bills
nopl <- filter(mdat, major!=21)
m.nopl <- glm(pv_all ~ majSD * minSD * any_statement + 
                majority + leader + les + mref + absPVI + 
                closerMinCS + elec.year + relSen + medDiff + majPriorityAll + 
                numSenate + sizeMajH, 
              family=binomial, data=filter(nopl, divided==1))

# R2s
r2s.pl <- sapply(list(m.full, m.nopl), 
                 function(x){
                   round(NagelkerkeR2(x)$R2, 3)
                 })

# Table A3: Regression table of full model and model without public lands bills
stargazer(m.full, m.nopl, type='latex', out='tableA3.tex',
          column.sep.width='0pt', no.space=TRUE, font.size='tiny', dep.var.labels.include = FALSE, dep.var.caption='', 
          title='Logistic regression models predicting whether a bill receives a final passage vote in the House.', model.numbers=FALSE,
          column.labels = c('Full Model', 'No Public Lands'),
          covariate.labels = c( # Main IV labels
            'Majority Spread', 'Minority Spread', 'Minority Priority',
            # Sponsor covariate labels
            'Majority', 'Leadership', 'Legislative Effectiveness', 'Member Referral Committee', 'Electoral Safety', 
            # Bill covariate labels
            'Closer to Minority', 'Election Year', 'Related in Senate', 'Party Difference', 'Majority Priority',
            # Political environment covariate labels
            'Seats in Senate', 'House Majority Size',
            # Interaction labels
            'Majority Spread * Minority Spread', 'Majority Spread * SOTU Mention', 'Minority Spread * SOTU Mention', 
            'Majority Spread * Minority Spread * SOTU Mention'
          ),
          keep.stat='n', add.lines = list(c('Nagelkerke Pseudo R2', r2s.pl))
)


### Figure A9: Predicted probabilities for full model and model without public lands bills

# Get data
pFull <- predDat(model=m.full, vars=c('any_statement', 'minSD', 'majSD'))
pNoPL <- predDat(model=m.nopl, vars=c('any_statement', 'minSD', 'majSD'))

#Plot
pdf("figureA9.pdf", width=10, height=10)
par(mfrow=c(2,2))
predProbPlot(pDat = pFull, breaks=3, ylab = 'Probability of Vote', xlab = 'Minority Spread', ylim=c(0.01, 0.08), var='SD',
             mains=c('Not a Minority Priority', 'Minority Priority'), rowControl = NULL, legTitle = 'Majority Spread')
predProbPlot(pDat = pNoPL, breaks = 3, ylab = 'Probability of Vote', xlab = 'Minority Spread', ylim=c(0.01, 0.08), var='SD',
             mains=c('Not a Minority Priority (No PL)', 'Minority Priority (No PL)'), rowControl = NULL, legTitle = 'Majority Spread')
par(mfrow=c(1,1))
dev.off()



###############################################################################################################################################################
### Figure A10: Changes in minority party priorities over time

### Format data
# Load, filter to years matching the congresses studied (1985-2006, though the 109th technically ends in 2007 it was before the 2007 sotu)
sotu <- fread(file='SOTU-topicmentions-1946-2019.csv', data.table=FALSE, stringsAsFactors=FALSE) %>%
  filter(year %in% 1985:2006) %>%
  dplyr::rename(minor=subtopic) 

# Match years to congresses
sotu$congress <- NA
congs <- 99:109
years <- sort(unique(sotu$year))
for(ii in 1:d(years)){
  if(years[ii] %% 2 == 1){ # If it's an odd year
    sotu$congress[sotu$year==years[ii]] <- congs[(ii+1)/2]
  } else {
    sotu$congress[sotu$year==years[ii]] <- congs[(ii/2)]
  }
}

# Drop years
sotu %<>% dplyr::select(-year)

# Create a Congress-by-subtopic matrix (one-hot encoded)
sotuMat <- matrix(data=0, nrow=d(unique(sotu$minor)), ncol=d(unique(sotu$congress)),
                  dimnames=list(paste('topic', unique(sotu$minor), sep='_'), 
                                paste('congress', unique(sotu$congress), sep='_'))
)
sm <- expand.grid(topic=sotu$minor, congress=sotu$congress)
congs <- unique(sm$congress)
for(ii in 1:d(congs)){
  df <- filter(sotu, congress==congs[ii])
  topics <- paste('topic', df$minor, sep='_')
  sotuMat[topics,ii] <- 1
}

# Shuffle the matrix by row (won't look exactly like in the paper)
sotuMat <- sotuMat[sample(nrow(sotuMat)),]

# Remove unified Congresses
unified <- paste(c(103, 107, 108, 109), collapse='|')
sotuMat <- sotuMat[,!str_detect(colnames(sotuMat), unified)]

# Plot in heatmap
pdf('figureA10.pdf', width=10, height=12)
heatmap(sotuMat, scale='none', col=c('black', 'white'), Rowv=NA, Colv=NA, xlab='Congresses with Divided Government', ylab='Topic',
        labRow=str_extract(rownames(sotuMat), '[0-9]+'),
        labCol=str_extract(colnames(sotuMat), '[0-9]+')
)
dev.off()



#########################################################################################################################################################################
### Figure A11: Using kurtosis instead of SD as a measure of spread

### Run models
# Vector of covariates for formulae
covs <- ' + majority + leader + les + mref + absPVI +
                closerMinCS + elec.year + relSen + medDiff + majPriorityAll +
                numSenate + sizeMajH'

# Model including SOTU mentions--adds motivation analysis 
m.sotu.kurt <- glm(paste0('pv_all ~ majKurt * minKurt * any_statement', covs), 
                   family=binomial, data=filter(mdat, divided==1))

# Models computing whether a bill became law (among all bills and those that passed the House)
m.law.kurt <- glm(paste0('plaw ~ majKurt * minKurt * any_statement', covs),
                  family=binomial, data=filter(mdat, divided==1))
m.law.pass.kurt <- glm(paste0('plaw ~ majKurt * minKurt * any_statement', covs),
                       family=binomial, data=filter(mdat, divided==1 & passh==1))

### Figure A11
# Get predicted data
pSotuKurt <- predDat(model=m.sotu.kurt, vars=c('any_statement', 'minKurt', 'majKurt'))
pLawKurt <- predDat(model=m.law.kurt, vars=c('any_statement', 'minKurt', 'majKurt'))
pLawPassKurt <- predDat(model=m.law.pass.kurt, vars=c('any_statement', 'minKurt', 'majKurt'))

# Pr(vote)
pdf("figureA11_vote.pdf", width=10, height=5)
predProbPlot(pDat = pSotuKurt, breaks = 3, ylab = 'Probability of Vote', xlab = 'Minority Kurtosis', ylim=c(0.01, 0.08), 
             mains=c('Not a Minority Priority', 'Minority Priority'), rowControl = c(1,2), legTitle = 'Majority Kurtosis',
             var='Kurt')
dev.off()

# Pr(law)
pdf("figureA11_law.pdf", width=10, height=5)
predProbPlot(pDat = pLawKurt, breaks = 3, ylab = 'Probability of Becoming Law', xlab = 'Minority Kurtosis', ylim=c(0, 0.05), 
             mains=c('Not a Minority Priority', 'Minority Priority'), rowControl = c(1,2), legTitle = 'Majority Kurtosis',
             var='Kurt')
dev.off()

# Pr(law|passed house)
pdf("figureA11_law_passed.pdf", width=10, height=5)
predProbPlot(pDat = pLawPassKurt, breaks = 3, ylab = 'Probability of Becoming Law', xlab = 'Minority Kurtosis', ylim=c(0.30, 0.70), 
             mains=c('Not a Minority Priority', 'Minority Priority'), rowControl = c(1,2), legTitle = 'Majority Kurtosis',
             var='Kurt')
dev.off()



################################################################################################################################################################
### Figures A12-A14: "Majority Priority" intead of "Majority Spread"

### Run models

# Full model (add interaction)
m.full.pri <- glm(pv_all ~ majPriorityAll * minSD + 
                majority + leader + les + mref + absPVI + 
                closerMinCS + elec.year + relSen + medDiff + majSD + 
                sizeMajH + numSenate + unityPty + presControl, 
              family=binomial, data=mdat)

# Model including SOTU mentions--adds motivation analysis 
# (only in cases of divided government between House and White House, which cuts out 103rd, 107th-109th congresses)
m.sotu.pri <- glm(pv_all ~ majPriorityAll * minSD * any_statement + 
                majority + leader + les + mref + absPVI + 
                closerMinCS + elec.year + relSen + medDiff + majSD +
                sizeMajH + numSenate + unityPty, 
              family=binomial, 
              data=filter(mdat, divided==1))

# Models computing whether a bill became law (among all bills and those that passed the House)
covs <- ' + majority + leader + les + mref + absPVI +
closerMinCS + elec.year + relSen + medDiff + majSD +
sizeMajH + numSenate + unityPty + presControl'
m.law.pri <- glm(paste0('plaw ~ majPriorityAll * minSD * any_statement', covs),
             family=binomial, data=filter(mdat, divided==1))
m.law.pass.pri <- glm(paste0('plaw ~ majPriorityAll * minSD * any_statement', covs),
                  family=binomial, data=filter(mdat, divided==1 & passh==1))


### Figure A12: 2-way interaction, Pr(vote)
# Get data
pDat <- predDat(model=m.full.pri, vars=c('minSD [all]', 'majPriorityAll'), data=mdat, std=TRUE)

# Plot
pdf("figureA12.pdf", width=5, height=5)

breaks <- 3
plot(0:(nrow(pDat)+2*breaks-1), 0:(nrow(pDat)+2*breaks-1), type='n', 
     ylab='Probability of Vote', xlab='Minority Spread', 
     ylim=c(0.02,0.05), xaxt='n', las=2, main='')
for(ii in 1:nrow(pDat)){
  colIndex <- which(c('No', 'Yes')==pDat$majPriorityAll[ii])
  if(ii %in% 1:2){
    points(ii, pDat$predicted[ii], col=cols3[colIndex+1], pch=16, cex=2)
    lines(x=rep(ii, 2), y=c(pDat$lo[ii], pDat$hi[ii]), col=cols3[colIndex+1], lwd=3)
  } else if(ii %in% 3:4){
    points(ii+2, pDat$predicted[ii], col=cols3[colIndex+1], pch=16, cex=2)
    lines(x=rep(ii+2, 2), y=c(pDat$lo[ii], pDat$hi[ii]), col=cols3[colIndex+1], lwd=3)
  } else if(ii %in% 5:6){
    points(ii+4, pDat$predicted[ii], col=cols3[colIndex+1], pch=16, cex=2)
    lines(x=rep(ii+4, 2), y=c(pDat$lo[ii], pDat$hi[ii]), col=cols3[colIndex+1], lwd=3)
  } 
}
axis(side=1, at=c(1.5, 5.5, 9.5), labels=c('Low', 'Median', 'High'), las=1, tick=FALSE)
legend(x='topright', legend=c('No', 'Yes'), fill=cols3[2:3], bty='n', title='Majority Priority')

dev.off()


## Figure A13: 3-way interaction, Pr(Vote)
p3 <- predDat(model=m.sotu.pri, vars=c('any_statement', 'minSD', 'majPriorityAll'))

# Set up plot parameters
pdf("figureA13.pdf", width=10, height=5)
par(mfrow=c(1,2))

p3No <- filter(p3, x==1)
breaks <- 3
plot(0:(nrow(pDat)+2*breaks-1), 0:(nrow(pDat)+2*breaks-1), type='n', 
     ylab='Probability of Vote', xlab='Minority Spread', 
     ylim=c(0.02,0.08), xaxt='n', las=2, main='Not a Minority Priority')
for(ii in 1:nrow(p3No)){
  colIndex <- which(c(0, 1)==p3No$majPriorityAll[ii])+1
  if(ii %in% 1:2){
    points(ii, p3No$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii, 2), y=c(p3No$lo[ii], p3No$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 3:4){
    points(ii+2, p3No$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+2, 2), y=c(p3No$lo[ii], p3No$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 5:6) {
    points(ii+4, p3No$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+4, 2), y=c(p3No$lo[ii], p3No$hi[ii]), col=cols3[colIndex], lwd=3)
  }
}
axis(side=1, at=c(1.5, 5.5, 9.5), labels=c('Low', 'Median', 'High'), las=1, tick=FALSE)
legend(x='topright', legend=c('No', 'Yes'), fill=cols3[2:3], bty='n', title='Majority Priority')

p3Yes <- filter(p3, x==2)
plot(0:(nrow(pDat)+2*breaks-1), 0:(nrow(pDat)+2*breaks-1), type='n', 
     ylab='Probability of Vote', xlab='Minority Spread', 
     ylim=c(0.02,0.08), xaxt='n', las=2, main='Minority Priority')
for(ii in 1:nrow(p3Yes)){
  colIndex <- which(c(0, 1)==p3Yes$majPriorityAll[ii])+1
  if(ii %in% 1:2){
    points(ii, p3Yes$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii, 2), y=c(p3Yes$lo[ii], p3Yes$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 3:4){
    points(ii+2, p3Yes$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+2, 2), y=c(p3Yes$lo[ii], p3Yes$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 5:6) {
    points(ii+4, p3Yes$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+4, 2), y=c(p3Yes$lo[ii], p3Yes$hi[ii]), col=cols3[colIndex], lwd=3)
  }
}
axis(side=1, at=c(1.5, 5.5, 9.5), labels=c('Low', 'Median', 'High'), las=1, tick=FALSE)
legend(x='topright', legend=c('No', 'Yes'), fill=cols3[2:3], bty='n', title='Majority Priority')

# Turn off graphics device
par(mfrow=c(1,1))
dev.off()


### Figure A14: 3-way interaction Pr(law)

# Get data
plaw <- predDat(model=m.law.pri, vars=c('any_statement', 'minSD', 'majPriorityAll'))

# Set up plot parameters
pdf("figureA14.pdf", width=10, height=5)
par(mfrow=c(1,2))

plawNo <- filter(plaw, x==1)
breaks <- 3
plot(0:(nrow(plawNo)+2*breaks-1), 0:(nrow(plawNo)+2*breaks-1), type='n', 
     ylab='Probability of Becoming Law', xlab='Minority Spread', 
     ylim=c(0,0.05), xaxt='n', las=2, main='Not a Minority Priority')
for(ii in 1:nrow(plawNo)){
  colIndex <- which(c(0, 1)==plawNo$majPriorityAll[ii])+1
  if(ii %in% 1:2){
    points(ii, plawNo$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii, 2), y=c(plawNo$lo[ii], plawNo$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 3:4){
    points(ii+2, plawNo$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+2, 2), y=c(plawNo$lo[ii], plawNo$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 5:6) {
    points(ii+4, plawNo$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+4, 2), y=c(plawNo$lo[ii], plawNo$hi[ii]), col=cols3[colIndex], lwd=3)
  }
}
axis(side=1, at=c(1.5, 5.5, 9.5), labels=c('Low', 'Median', 'High'), las=1, tick=FALSE)
legend(x='topright', legend=c('No', 'Yes'), fill=cols3[2:3], bty='n', title='Majority Priority')

plawYes <- filter(plaw, x==2)
plot(0:(nrow(plawYes)+2*breaks-1), 0:(nrow(plawYes)+2*breaks-1), type='n', 
     ylab='Probability of Becoming Law', xlab='Minority Spread', 
     ylim=c(0, 0.05), xaxt='n', las=2, main='Minority Priority')
for(ii in 1:nrow(plawYes)){
  colIndex <- which(c(0, 1)==plawYes$majPriorityAll[ii])+1
  if(ii %in% 1:2){
    points(ii, plawYes$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii, 2), y=c(plawYes$lo[ii], plawYes$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 3:4){
    points(ii+2, plawYes$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+2, 2), y=c(plawYes$lo[ii], plawYes$hi[ii]), col=cols3[colIndex], lwd=3)
  } else if(ii %in% 5:6) {
    points(ii+4, plawYes$predicted[ii], col=cols3[colIndex], pch=16, cex=2)
    lines(x=rep(ii+4, 2), y=c(plawYes$lo[ii], plawYes$hi[ii]), col=cols3[colIndex], lwd=3)
  }
}
axis(side=1, at=c(1.5, 5.5, 9.5), labels=c('Low', 'Median', 'High'), las=1, tick=FALSE)
legend(x='topright', legend=c('No', 'Yes'), fill=cols3[2:3], bty='n', title='Majority Priority')

# Turn off graphics device
par(mfrow=c(1,1))
dev.off()
